13.2 GAS与输入

前文中在蓝图中使用输入事件配合TryActivateAbility来实现技能的触发是一种比较初级的技能调用方案,目前作为触发的方案来讲并不会有太大的问题。但其实质上并没有在内部将输入与技能绑定起来,具体来说,GA中有一些监听输入的任务,如WaitForInputRelease,并不会响应我们目前简单绑定的输入响应。

在这里插入图片描述
在这里插入图片描述

如图,我们对GA_FirstAbility做出简单修改,添加了一个输入相关的异步任务。从打印结果,可见WaitInputRelease的异步任务并没有正常响应。我们希望在能够直接将输入系统中设定的输入和对应的技能进行绑定。同时在调试调试信息中显示,创建的技能和异步任务也一直存在于内存中没有正确释放:

在这里插入图片描述
简单来说,本文的最终结果就是令GA在包含上面的逻辑功能时能够正常执行并且结束。

在UE的最新版本中引入了Enhanced Input(基于原输入系统拓展的一套进阶输入控制方案),本文将分别实现两种输入系统的输入 - GA的绑定。

1. 输入测试环境搭建

先将原先创建的通过T来TryActivateAbility的蓝图内容删除(以及GiveAbility,后续会给出更合适的写法和位置),防止干扰我们的后续内容。

针对旧输入系统的配置:

  1. 在输入配置中添加按键绑定:OldTestInput - T;
  2. 在代码中添加枚举EGASAbilityInputID
  3. 创建专为旧输入系统测试的GA:GA_OldInputBindAbility,直接继承前面修改过的GA_FirstAbility即可,无须其他逻辑;

在这里插入图片描述

UENUM(BlueprintType)
enum class EGASAbilityInputID : uint8
{
	// 0 None
	None			UMETA(DisplayName = "None"),

	// 1 Jump
	Jump			UMETA(DisplayName = "Jump"),
	
	// 2 T
	OldTestInput		UMETA(DisplayName = "OldTestInput"),
};

在这里插入图片描述

针对EnhancedInput输入系统的配置:

  1. 在插件界面,启用EnhancedInput,并且重启编辑器;
  2. 在输入配置界面,将输入相关的默认类切换为EnhancedInput的内容(EnhancedInput的类是继承自默认输入(旧输入)系统的,所以在启用EnhancedInput时,不会干扰旧的输入,及旧输入系统定义的按键绑定也可以生效);
  3. 创建专为EnhancedInput输入系统测试的InputAction:IA_EnhancedInputTest,并将其添加到创建的默认的InputMappingContext里,设定按键映射到E,即IMC_Default - IA_EnhancedInputTest - E;
  4. 创建专为EnhancedInput输入系统测试的GA:GA_EnhancedInputBindAbility,直接继承前面修改过的GA_FirstAbility即可,无须其他逻辑;
  5. 项目名.build.cs文件中添加EnhancedInput的模块依赖;

在这里插入图片描述
在这里插入图片描述

PrivateDependencyModuleNames.AddRange(new string[]
{
	"GameplayAbilities",
	"GameplayTasks",
	"GameplayTags",
	"EnhancedInput"
});

在这里插入图片描述

2. GA和旧输入系统的绑定方案

在角色类GASCharacterBase的头文件中,添加技能和输入的容器:

UPROPERTY(EditDefaultsOnly, BlueprintReadOnly, Category = "GAS|Binding")
TMap<EGASAbilityInputID, TSubclassOf<UGameplayAbility>> OldInputDefaultAbilities;

这里就是负责保存输入和技能映射的位置。

在该类PossessedBy方法中,在初始化技能系统组件之后,加入技能的赋予相关代码:

if (IsValid(AbilitySystemComponent))
{
	AbilitySystemComponent->InitAbilityActorInfo(this, this);

	if (GetLocalRole() == ROLE_Authority)
	{
		for (TTuple<EGASAbilityInputID, TSubclassOf<UGameplayAbility>>& Ability : OldInputDefaultAbilities)
		{
			AbilitySystemComponent->GiveAbility(
				FGameplayAbilitySpec(Ability.Value, 1, static_cast<int32>(Ability.Key), this));
		}
	}
}

这样,在游戏开始后(具体是角色被控制器控制时),技能就被添加给角色(角色身上的技能系统组件)。

最后进行技能和输入之间的绑定,推荐的写的位置是SetupPlayerInputComponent方法中(最好是在GASCharacterBase的派生类中),这个方法会在前面提到的PossessedBy之后调用,确保技能系统组件已经正确初始化并且技能已经被赋予。

if (IsValid(AbilitySystemComponent))
{
	AbilitySystemComponent->BindAbilityActivationToInputComponent(InputComponent,
	                                                              FGameplayAbilityInputBinds(
		                                                              FString(), FString(),
		                                                              FString("EGASAbilityInputID"), -1, -1));
}

注意: 这里的EGASAbilityInputID就是之前定义好的输入的枚举,这个函数内部是依赖于字符串匹配的,所以需要注意拼写。

编译运行,进入到角色蓝图中,将已经预制好的输入和GA填入:

在这里插入图片描述
结果:
按下T键
在这里插入图片描述
松开T键
在这里插入图片描述
技能及技能任务都正确触发。

这里整套的旧输入与技能的绑定方案参考并简化了GASDocumentation中的实现。

3. GA和新输入系统(Enhanced Input)的绑定方案

旧的输入绑定方案是建立在GAS本就对其进行了一定程度的支持的基础(主要是BindAbilityActivationToInputComponent)之上,所以在接入EnhancedInput后,我们需要稍稍对AbilitySystemComponent进行一些扩展,在实现功能的前提下,对函数接口进行统一(即新方案可以兼容旧方案,两者可以同时使用且不会有冲突)。

创建EnhancedAbilitySystemComponent类,继承AbilitySystemComponent,后面也将使用自定义扩展的技能系统组件:

#pragma once

#include "CoreMinimal.h"
#include "AbilitySystemComponent.h"
#include "InputAction.h"
#include "EnhancedAbilitySystemComponent.generated.h"

USTRUCT()
struct FAbilityInputBinding
{
	GENERATED_BODY()

	int32 InputID = 0;
	uint32 OnPressedHandle = 0;
	uint32 OnReleasedHandle = 0;
	TArray<FGameplayAbilitySpecHandle> BoundAbilitiesStack;
};

UCLASS()
class INSIDEGAS_API UEnhancedAbilitySystemComponent : public UAbilitySystemComponent
{
	GENERATED_BODY()

public:
	// Sets default values for this component's properties
	UEnhancedAbilitySystemComponent(const FObjectInitializer& ObjectInitializer);

	virtual void BindAbilityActivationToInputComponent(UInputComponent* InputComponent,
	                                                   FGameplayAbilityInputBinds BindInfo) override;

	void BindAbilityActivationToInputComponent(UEnhancedInputComponent* EnhancedInputComponent,
	                                           TMap<UInputAction*, TSubclassOf<UGameplayAbility>> InputAbilityMap);

	UFUNCTION()
	void OnAbilityBindedInputPressed(UInputAction* InputAction);

	UFUNCTION()
	void OnAbilityBindedInputReleased(UInputAction* InputAction);

	UPROPERTY(Transient)
	TMap<UInputAction*, FAbilityInputBinding> MappedAbilities;
};
#include "EnhancedAbilitySystemComponent.h"
#include "EnhancedInputComponent.h"

namespace AbilityInputBindingComponent_Impl
{
	constexpr int32 InvalidInputID = 10;
	int32 IncrementingInputID = InvalidInputID;

	static int32 GetNextInputID()
	{
		return ++IncrementingInputID;
	}
}

UEnhancedAbilitySystemComponent::UEnhancedAbilitySystemComponent(const FObjectInitializer& ObjectInitializer)
	: Super(ObjectInitializer)
{
}

void UEnhancedAbilitySystemComponent::BindAbilityActivationToInputComponent(UInputComponent* InputComponent,
	FGameplayAbilityInputBinds BindInfo)
{
	Super::BindAbilityActivationToInputComponent(InputComponent, BindInfo);
}

void UEnhancedAbilitySystemComponent::BindAbilityActivationToInputComponent(
	UEnhancedInputComponent* EnhancedInputComponent, TMap<UInputAction*, TSubclassOf<UGameplayAbility>> InputAbilityMap)
{
	if (ensureMsgf(EnhancedInputComponent,
	               TEXT("Project must use EnhancedInputComponent to support PlayerControlsComponent")))
	{
		for (auto& InputBinding : MappedAbilities)
		{
			const int32 NewInputID = AbilityInputBindingComponent_Impl::GetNextInputID();
			InputBinding.Value.InputID = NewInputID;

			for (FGameplayAbilitySpecHandle AbilityHandle : InputBinding.Value.BoundAbilitiesStack)
			{
				FGameplayAbilitySpec* FoundAbility = FindAbilitySpecFromHandle(AbilityHandle);
				if (FoundAbility != nullptr)
				{
					FoundAbility->InputID = NewInputID;
				}
			}
		}
		
		for (auto InputAbility : InputAbilityMap)
		{
			FGameplayAbilitySpec AbilitySpec = BuildAbilitySpecFromClass(InputAbility.Value);

			if (IsValid(AbilitySpec.Ability))
			{
				FGameplayAbilitySpecHandle AbilitySpecHandle = GiveAbility(AbilitySpec);
				
				using namespace AbilityInputBindingComponent_Impl;

				FGameplayAbilitySpec* BindingAbility = FindAbilitySpecFromHandle(AbilitySpecHandle);

				UInputAction* InputAction = InputAbility.Key;

				FAbilityInputBinding* AbilityInputBinding = MappedAbilities.Find(InputAction);
				
				if (AbilityInputBinding)
				{
					FGameplayAbilitySpec* OldBoundAbility = FindAbilitySpecFromHandle(AbilityInputBinding->BoundAbilitiesStack.Top());
					if (OldBoundAbility && OldBoundAbility->InputID == AbilityInputBinding->InputID)
					{
						OldBoundAbility->InputID = InvalidInputID;
					}
				}
				else
				{
					AbilityInputBinding = &MappedAbilities.Add(InputAction);
					AbilityInputBinding->InputID = GetNextInputID();
				}

				if (BindingAbility)
				{
					BindingAbility->InputID = AbilityInputBinding->InputID;
				}

				AbilityInputBinding->BoundAbilitiesStack.Push(AbilitySpecHandle);

				// Pressed event
				if (AbilityInputBinding->OnPressedHandle == 0)
				{
					AbilityInputBinding->OnPressedHandle = EnhancedInputComponent->BindAction(InputAction, ETriggerEvent::Started, this, &UEnhancedAbilitySystemComponent::OnAbilityBindedInputPressed, InputAction).GetHandle();
				}

				// Released event
				if (AbilityInputBinding->OnReleasedHandle == 0)
				{
					AbilityInputBinding->OnReleasedHandle = EnhancedInputComponent->BindAction(InputAction, ETriggerEvent::Completed, this, &UEnhancedAbilitySystemComponent::OnAbilityBindedInputReleased, InputAction).GetHandle();
				}
			}
		}
	}
}

void UEnhancedAbilitySystemComponent::OnAbilityBindedInputPressed(UInputAction* InputAction)
{
	using namespace AbilityInputBindingComponent_Impl;

	FAbilityInputBinding* FoundBinding = MappedAbilities.Find(InputAction);
	if (FoundBinding && ensure(FoundBinding->InputID != InvalidInputID))
	{
		AbilityLocalInputPressed(FoundBinding->InputID);
	}
}

void UEnhancedAbilitySystemComponent::OnAbilityBindedInputReleased(UInputAction* InputAction)
{
	using namespace AbilityInputBindingComponent_Impl;

	FAbilityInputBinding* FoundBinding = MappedAbilities.Find(InputAction);
	if (FoundBinding && ensure(FoundBinding->InputID != InvalidInputID))
	{
		AbilityLocalInputReleased(FoundBinding->InputID);
	}
}

其中BindAbilityActivationToInputComponent是核心函数接口,也是主要进行重载的函数,使其可以对扩展输入EnhancedInputComponent进行绑定。

注意: 另外InvalidInputID是我们所使用的InputID的起始位置,如果设定其为0的话,会跟旧输入的绑定相冲突,这里建议设定一定的偏移(我这里设定为10,可以根据使用的旧输入的数量进行适当的修改)。

到角色部分进行绑定的调用,不过首先要进行默认组件的修改:

AInsideGASCharacter::AInsideGASCharacter(const FObjectInitializer& ObjectInitializer)
	: Super(ObjectInitializer.SetDefaultSubobjectClass<UEnhancedAbilitySystemComponent>(
		AGASCharacterBase::AbilitySystemComponentName))

然后到SetupInputComponent中添加输入绑定代码:

EnhancedAbilitySystemComponent = Cast<UEnhancedAbilitySystemComponent>(AbilitySystemComponent);
UEnhancedInputComponent* EnhancedInputComponent = Cast<UEnhancedInputComponent>(InputComponent);

UEnhancedInputLocalPlayerSubsystem* Subsystem = GetController<APlayerController>()->GetLocalPlayer()->GetSubsystem<UEnhancedInputLocalPlayerSubsystem>();
check(Subsystem);

if (DefaultInputMappingContext)
{
	Subsystem->AddMappingContext(DefaultInputMappingContext, InputPriority);
}

if (IsValid(EnhancedAbilitySystemComponent) && IsValid(EnhancedInputComponent))
{
	EnhancedAbilitySystemComponent->BindAbilityActivationToInputComponent(
		EnhancedInputComponent, EnhancedInputDefaultAbilities);
}

进入到角色蓝图中,挂载引用:

在这里插入图片描述
结果:
按下E键

在这里插入图片描述
松开E键

在这里插入图片描述
功能需求满足。

这里的方案参考了古代山谷项目中的实现方法,古代山谷和Lyra项目中使用GameFeature(具体是FeatureAction的AddAbilities)进行调用,并将功能整合在一个额外的PawnComponent上。我在本文中对这个方案进行了一定程度的简化,摒弃了GameFeature的调用,而直接在角色身上进行绑定。好处是简单直观,但是相比之下,缺少了一定的灵活性。

4. 总结

本文从应用层面对GAS的输入触发进行了一定程度的归纳,并分别针对旧/新输入系统提出了与GAS结合使用的技术方案。需要注意的是,这里的方案不建议直接应用到实际项目,其中还缺少了一些生命周期相关的逻辑(本文关注激活与绑定,并没有处理卸载和解绑)。

后续的文章将采用最新的EnhancedInput输入系统,输入的绑定方案也将采用上述的第二种。

参考:

虚幻5新特性之EnhancedInput

游戏技能 - Gameplay Abilities

请在场景1代码的基础上修改场景2,优化变量全部改成场景1的,请给出修改后的场景2代码,即场景1补充CCS后的代码:%% 场景1:考虑碳交易+P2G clc clear %% 加载电负荷、热负荷 load date.mat; P_load_e = date(1, :); % 总电负荷 H_h = date(2, :); % 高品位热负荷 H_l = date(3, :); % 低品位热负荷 H_load_h = H_h + H_l; % 总热负荷 % 成本系数和能源价格 c_coal = [0.0004818, 2.695, 52, 0.0004751, 0.3536, 0.001004] * 0.1471; p_coal = 800; % 标煤价格 (元/吨) p_gas = 3.2; % 天然气价格 (元/m³) p_buy_e = 1000 * [0.17, 0.17, 0.17, 0.17, 0.17, 0.17, 0.17, 0.49, 0.83, 0.83, 0.83, 0.83, 0.49, 0.49, 0.49, 0.49, 0.49, 0.83, 0.83, 0.83, 0.83, 0.49, 0.49, 0.49]; % 电价系数 V_buy_h = 5 * [10, 10.4, 11.3, 11.7, 12.4, 13.2, 14.1, 14.6, 15.1, 15.3, 14.7, 14.4, 14, 13.8, 13.3, 12.7, 12.2, 12.3, 11.9, 11.8, 11.4, 11.1, 10.4, 10.2]; % 热价基数 C_sell_h = 180; % 售热价格 (元/MWh) % 气象数据(辐照度 kW/m&sup2;,风速 m/s) P_Solar = [0,0,0,0,0,0.09,1.15,2.82,4.54,6.20,7.44,8.01,7.93,7.09,5.93,4.57,2.97,1.27,0.15,0,0,0,0,0]; P_Wind = [7.14,7.20,7.22,6.89,7.18,7.07,6.43,5.97,6.33,6.56,6.63,6.70,6.65,6.92,6.86,6.85,6.89,7.01,7.02,6.97,7.05,7.06,7.03,7.10]; %% 新能源模型参数 % 光伏模型 (单位: MW) S_PV = 40000 * 3.75; % 光伏面积 (m&sup2;) k_eta_PV = 0.157; % 光电转换效率 P_PV_max = zeros(1, 24); for i = 1:24 P_PV_max(i) = (S_PV * P_Solar(i) * k_eta_PV) / 1000; end % 风电模型 (单位: MW) Q_S_WT = 80000 * 4.8; % 单机额定容量 (kW) v_cut_in = 3; v_rated = 12; v_cut_out = 20; % 风速参数 P_WT_max = zeros(1, 24); for i = 1:24 v = P_Wind(i); if v >= v_cut_in && v <= v_rated P_WT_max(i) = Q_S_WT * (v - v_cut_in) / (v_rated - v_cut_in) / 1000; elseif v > v_rated && v <= v_cut_out P_WT_max(i) = Q_S_WT / 1000; else P_WT_max(i) = 0; end end %% 定义决策变量 P_th_h = sdpvar(1, 24, 'full'); % 火电机组热出力 P_th_e = sdpvar(1, 24, 'full'); % 火电机组电出力 z_th_1 = binvar(1, 24); % 机组1状态 z_th_2 = binvar(1, 24); % 机组2状态 P_e_eb = sdpvar(1, 24, 'full'); % 电锅炉电功率 P_g_gb = sdpvar(1, 24, 'full'); % 燃气锅炉燃气功率 P_g_chp = sdpvar(1, 24, 'full'); % CHP机组燃气功率 B_cr_h = binvar(1, 24); % 储热状态 (0:储热, 1:放热) P_loss_e = sdpvar(1, 24, 'full'); % 失负荷功率 m_ccs_st = sdpvar(1, 24, 'full'); % CCS解吸量 m_ccs_ab = sdpvar(1, 24, 'full'); % CCS吸收量 P_e_p2g = sdpvar(1, 24, 'full'); % P2G用电功率 P_h2_chp = sdpvar(1, 24, 'full'); % CHP掺氢功率 level_cool = binvar(3, 24); % COAP冷却级别 Z_coap_mt = sdpvar(1, 24, 'full'); % 活性炭质量 P_e_coap = sdpvar(1, 24, 'full'); % COAP能耗 B_coap_x1 = sdpvar(2, 24, 'full'); % COAP辅助变量 P_pv_e = sdpvar(1, 24, 'full'); % 光伏上网功率 P_wt_e = sdpvar(1, 24, 'full'); % 风电上网功率 P_ess_e = sdpvar(1, 24, 'full'); % 电化学储能净功率 P_ess_cha = sdpvar(1, 24, 'full'); % 储能充电功率 P_ess_dis = sdpvar(1, 24, 'full'); % 储能放电功率 H_tes_dis = sdpvar(1, 24, 'full'); % 储热罐释热功率 H_tes_cha = sdpvar(1, 24, 'full'); % 储热罐储热功率 SOC_HSS = sdpvar(1,24); %% 系统参数 SOC_HSS_max = 80; % 储热罐最大容量 (MWh) SOC_ESS_max = 40; % 储电能最大容量 (MWh) P_eb_max = 10; % 电锅炉最大功率 (MW) P_gb_max = 40; % 燃气锅炉最大功率 (MW) k_eb = 0.99; % 电锅炉效率 k_gb = 0.81; % 燃气锅炉效率 K_chp_h = 0.47; % CHP机组气-热比 K_chp_e = 0.35; % CHP机组气-电比 %% 约束条件 n = 24; Cons = []; for i = 1:n % 火电热功率约束 Cons = [Cons, (0 <= P_th_h(1,i)) & (P_th_h(1,i) <= 2*98.18)]; Cons = [Cons, implies(P_th_h(1,i) >= 2*57.27, z_th_1(1,i) == 1)]; Cons = [Cons, z_th_1(1,i) + z_th_2(1,i) == 1]; % 火电电功率约束 Cons = [Cons, (2*43.37 <= P_th_e(1,i)) & (P_th_e(1,i) <= 2*135)]; % 火电热-电耦合约束 Cons = [Cons, implies(z_th_2(1,i), (P_th_e(1,i) <= 2*135 - 0.25*P_th_h(1,i)) & (P_th_e(1,i) >= 2*54 - 0.1857*P_th_h(1,i)))]; Cons = [Cons, implies(z_th_1(1,i), (P_th_e(1,i) <= 2*135 - 0.25*P_th_h(1,i)) & (P_th_e(1,i) >= 1.64*P_th_h(1,i) - 2*50.56))]; % CHP机组约束 Cons = [Cons, (0 <= P_g_chp(1,i)) & (P_g_chp(1,i) <= 2*50)]; % 电锅炉约束 Cons = [Cons, (0 <= P_e_eb(1,i)) & (P_e_eb(1,i) <= P_eb_max)]; % 失负荷约束 Cons = [Cons, P_loss_e(1,i) >= 0]; % 储热系统功率约束 Cons = [Cons, (0 <= H_tes_cha(1,i)) & (H_tes_cha(1,i) <= 20)]; Cons = [Cons, (0 <= H_tes_dis(1,i)) & (H_tes_dis(1,i) <= 20)]; % 储热容量约束 (动态) if i == 1 initial_H = 0.3 * SOC_HSS_max; else initial_H = 0; end Cons = [Cons, (0 <= initial_H + sum(0.98*H_tes_cha(1,1:i) + k_eb*P_e_eb(1,1:i) - H_tes_dis(1,1:i)) <= SOC_HSS_max)]; % 储热状态互斥 Cons = [Cons, implies(B_cr_h(i) == 0, H_tes_dis(1,i) == 0)]; Cons = [Cons, implies(B_cr_h(i) == 1, H_tes_cha(1,i) == 0)]; % 电功率平衡 Cons = [Cons, (P_th_e(1,i) + K_chp_e*P_g_chp(1,i) + P_wt_e(1,i) + P_pv_e(1,i) + P_loss_e(1,i) - P_e_eb(1,i) - P_e_p2g(1,i) - P_e_coap(1,i) + P_ess_e(1,i) - (5 + 0.269*m_ccs_st(1,i))) == P_load_e(i)]; % 热功率平衡 Cons = [Cons, P_th_h(1,i) + 0.98*H_tes_dis(1,i) - H_tes_cha(1,i) + K_chp_h*P_g_chp(1,i) + k_gb*P_g_gb(1,i) == H_load_h(i)]; % 燃气锅炉约束 Cons = [Cons, (0 <= P_g_gb(1,i)) & (P_g_gb(1,i) <= P_gb_max)]; % CCS约束 (当前禁用) Cons = [Cons, m_ccs_st(1,i) == 0]; Cons = [Cons, m_ccs_ab(1,i) == 0]; % P2G约束 Cons = [Cons, (0 <= P_e_p2g(1,i)) & (P_e_p2g(1,i) <= 30)]; Cons = [Cons, (0 <= P_h2_chp(1,i)) & (P_h2_chp(1,i) <= P_e_p2g(1,i))]; Cons = [Cons, (P_g_chp(1,i) + P_g_gb(1,i))*0.003539 >= 5*P_h2_chp(1,i)*0.01083]; % 掺氢比约束 % COAP冷却级别 (固定使用第三级) Cons = [Cons, level_cool(1,i) + level_cool(2,i) + level_cool(3,i) == 1]; Cons = [Cons, level_cool(3,i) == 1]; Cons = [Cons, P_e_coap(1,i) == 0]; % 所有级别能耗均为0 % 活性炭约束 Cons = [Cons, (0 <= Z_coap_mt(1,i)) & (Z_coap_mt(1,i) <= 20)]; % 电化学储能约束 Cons = [Cons, (-20 <= P_ess_e(1,i)) & (P_ess_e(1,i) <= 20)]; Cons = [Cons, (-20 <= P_ess_cha(1,i)) & (P_ess_cha(1,i) <= 0)]; Cons = [Cons, (0 <= P_ess_dis(1,i)) & (P_ess_dis(1,i) <= 20)]; Cons = [Cons, implies(P_ess_e(1,i) >= 0, (P_ess_dis(1,i) == P_ess_e(1,i)) & (P_ess_cha(1,i) == 0))]; Cons = [Cons, implies(P_ess_e(1,i) <= 0, (P_ess_dis(1,i) == 0) & (P_ess_cha(1,i) == P_ess_e(1,i)))]; % 储能容量约束 (动态) if i == 1 initial_ESS_SOC = 0.3 * SOC_ESS_max; else initial_ESS_SOC = 0; end Cons = [Cons, (0 <= initial_ESS_SOC + sum(P_ess_dis(1,1:i))/0.95 + 0.95*sum(P_ess_cha(1,1:i)) <= SOC_ESS_max)]; % COAP辅助变量约束 Cons = [Cons, B_coap_x1(1,i) == (P_th_e(1,i) + P_th_h(1,i)) * 6.48]; Cons = [Cons, B_coap_x1(2,i) == (P_th_e(1,i) + P_th_h(1,i)) * 2.88]; % 新能源出力约束 Cons = [Cons, (0 <= P_pv_e(1,i)) & (P_pv_e(1,i) <= P_PV_max(i))]; Cons = [Cons, (0 <= P_wt_e(1,i)) & (P_wt_e(1,i) <= P_WT_max(i))]; end % 燃气轮机爬坡约束 for i = 2:n Cons = [Cons, (2*-28.6 <= P_g_chp(1,i) - P_g_chp(1,i-1)) & (P_g_chp(1,i) - P_g_chp(1,i-1) <= 2*28.6)]; end % 储热和储能周期平衡约束 Cons = [Cons, (-1 <= sum(0.98*H_tes_cha(1,1:24) + k_eb*P_e_eb(1,1:24) - H_tes_dis(1,1:24)) <= 1)]; Cons = [Cons, (-1 <= sum(P_ess_dis(1,1:24))/0.95 + 0.95*sum(P_ess_cha(1,1:24)) <= 1)]; %% 燃料购能成本 Fun = 0; for i = 1:n % 火电机组燃料成本 C_th_c = p_coal * (c_coal(3) + c_coal(2)*P_th_e(i) + c_coal(5)*P_th_h(i) + c_coal(1)*P_th_e(i)^2 + c_coal(6)*P_th_e(i)*P_th_h(i) + c_coal(4)*P_th_h(i)^2); % CHP机组和燃气锅炉燃气成本 C_chp_c = p_gas * (P_g_chp(i) / 0.01083); C_gb_c = p_gas * (P_g_gb(i) / 0.01083); % 环保代价成本 C_cb_c = 6.4 * B_coap_x1(1,i) + 8 * B_coap_x1(2,i) + 3.16; % 碳排放计算 C_mc = 2.57 * (C_th_c/p_coal) + 0.054*0.98*3.67*(P_g_chp(i) + P_g_gb(i) - P_h2_chp(i)); % 碳配额计算 C_Ec = 0.69135*(P_th_e(i) + K_chp_e*P_g_chp(i)) + 0.3*(P_th_h(i) + K_chp_h*P_g_chp(i) + k_gb*P_g_gb(i)); % 碳交易成本 C_ct_c = 205 * (C_mc - C_Ec); % 失负荷惩罚成本 C_buy_c = 5 * p_buy_e(i) * P_loss_e(i); % 售热售电收益 I_sale = p_buy_e(i) * (P_th_e(i) + K_chp_e*P_g_chp(i) + P_wt_e(i) + P_pv_e(i) - P_e_eb(i) - P_e_p2g(i) - P_e_coap(i) + P_ess_e(i)) + C_sell_h * (P_th_h(i) + 0.98*H_tes_dis(i) - H_tes_cha(i) + K_chp_h*P_g_chp(i) + k_gb*P_g_gb(i)); % P2G产气收益 I_ch4_c = p_gas * (((P_e_p2g(i)*0.85 - P_h2_chp(i)) * 0.75) / 0.01083); % 总成本 C0_all = C_th_c + C_gb_c + C_chp_c - p_gas*P_h2_chp(i)/(0.01083) + C_cb_c + C_buy_c + C_ct_c - I_ch4_c; % 目标函数累加 Fun = Fun + (I_sale - C0_all); end %% 设备投资和运维成本 % 火电建设成本 C_Ld = (4589000*((135*2)*(1+0.07)) + 100) / (365*30); C_cs = 0; for i=1:n % 累加燃料相关成本 fuel_cost_i = C_th_c + C_gb_c + C_chp_c - p_gas*P_h2_chp(i)/0.01083 - I_ch4_c; C_cs = C_cs + fuel_cost_i; end % 火电运维成本 C_com = (0.1536/(1-0.1536)) * (C_Ld + C_cs); % 光伏成本(建设+运维) C_pv = (6000*1000*200)/(365*25) + 12*sum(P_PV_max(1:24)); % 风电成本(建设+运维) C_wt = (8000*1000*440)/(365*20) + 20*sum(P_WT_max(1:24)); % 电锅炉成本 C_eb = (3200*1000*P_eb_max)/(365*20) + 2*sum(P_e_eb(1,1:24)); % 储热罐成本 C_yan = (300*1000*SOC_HSS_max)/(365*5) + 1.5*sum(H_tes_cha(1,1:24)); % 电化学储能成本 C_ess = (2000*1000*SOC_ESS_max)/(365*5); % 煤电容量电价补偿 I_c_e = (330*1000*(370))/365; % 总设备成本 C_q_om = C_Ld + C_com + C_pv + C_wt + C_eb + C_yan + C_ess - I_c_e; %% 弃风弃光惩罚 P_v_a = P_PV_max - P_pv_e; % 光伏弃光量 P_w_a = P_WT_max - P_wt_e; % 风电弃风量 F_a_em = 257.396 * sum(P_w_a + P_v_a); % 惩罚成本 %% 目标函数 Fun = Fun - C_q_om - F_a_em; % 总收益减去设备成本和弃能惩罚 Fun = -Fun; % 转换为最小化问题 %% CPLEX优化求解 ops = sdpsettings('solver', 'cplex', 'verbose', 2); ops.cplex.timelimit = 120; % 设置求解时间限制为120秒 optimize(Cons, Fun, ops); %% 结果可视化 figure(1) x = 1:24; hold on; P_pv_e_val = value(P_pv_e); P_wt_e_val = value(P_wt_e); % 初始化弃能区域标记 has_pv_curtail = false; has_wt_curtail = false; % 绘制弃能区域 for i = 1:24 % 光伏弃光区域 if P_PV_max(i) > P_pv_e_val(i) patch([x(i)-0.4, x(i)+0.4, x(i)+0.4, x(i)-0.4], ... [P_pv_e_val(i), P_pv_e_val(i), P_PV_max(i), P_PV_max(i)], ... [0.97, 0.75, 0.40], 'EdgeColor', 'none', 'FaceAlpha', 0.5); has_pv_curtail = true; end % 风电弃风区域 if P_WT_max(i) > P_wt_e_val(i) patch([x(i)-0.4, x(i)+0.4, x(i)+0.4, x(i)-0.4], ... [P_wt_e_val(i), P_wt_e_val(i), P_WT_max(i), P_WT_max(i)], ... [0.40, 0.70, 0.95], 'EdgeColor', 'none', 'FaceAlpha', 0.5); has_wt_curtail = true; end end % 绘制实际出力曲线 p_pv_actual = plot(x, P_pv_e_val, '-o', 'Color', [0.90, 0.55, 0.20], 'LineWidth', 1, 'MarkerSize', 5); p_wt_actual = plot(x, P_wt_e_val, '-o', 'Color', [0.20, 0.50, 0.80], 'LineWidth', 1, 'MarkerSize', 5); % 绘制最大出力曲线 p_pv_max = plot(x, P_PV_max, '-s', 'Color', [0.80, 0.35, 0.10], 'LineWidth', 1); p_wt_max = plot(x, P_WT_max, '-s', 'Color', [0.10, 0.30, 0.60], 'LineWidth', 1); % 创建图例项 legend_items = [p_pv_actual, p_pv_max]; legend_labels = {'Actual PV Output', 'Maximum PV Output'}; if has_pv_curtail legend_items(end+1) = patch(NaN, NaN, [0.97, 0.75, 0.40], 'FaceAlpha', 0.5); legend_labels{end+1} = 'Curtailed PV Power'; end legend_items(end+1:end+2) = [p_wt_actual, p_wt_max]; legend_labels(end+1:end+2) = {'Actual WT Output', 'Maximum WT Outpu'}; if has_wt_curtail legend_items(end+1) = patch(NaN, NaN, [0.40, 0.70, 0.95], 'FaceAlpha', 0.5); legend_labels{end+1} = 'Curtailed WT Power'; end legend(legend_items, legend_labels); grid on; xlabel('time (h)'); ylabel('Power (MW)'); title('Output Curve of New Energy Units'); xlim([0, 25]); ylim([-1 260]); set(gca, 'XTick', 0:5:25); box on; hold off; figure(2) % Electric Power Balance subplot(1,3,1); hold on; P = [P_g_chp*K_chp_e; P_th_e; P_pv_e; P_wt_e; -P_e_eb; -P_e_p2g; P_ess_cha; P_ess_dis]'; e = bar(P, 'stacked'); hold on; plot(1:24, P_load_e, '-g*'); grid on; colors_e = [ 0.7 0.2 0.7 % CHP (purple) 0.9 0.4 0.2 % Thermal (orange) 1 0.8 0.2 % PV (yellow) 0.5 0.8 0.2 % WT (green) 0.6 0.6 0.6 % EB (gray) 0.8 0.3 0.3 % P2G (red) 0.2 0.8 0.8 % ESS Charge (cyan) 0.3 0.7 0.9 % ESS Discharge (light blue) ]; hold on; for i = 1:size(colors_e,1) e(i).FaceColor = colors_e(i,:); e(i).EdgeColor = 'none'; end hold on; legend_items = {'CHP','CTP','PV','WT','EB','P2G','ESS(cha)','ESS(dis)','Electric Load'}; legend(legend_items); xlabel('Time (h)'); ylabel('Power (MW)'); axis([0,25,-200,700]); xticks(0:2:24); title('Electric Power Balance'); % Thermal Power Balance subplot(1,3,2); H1 = [K_chp_h*P_g_chp; P_th_h; k_gb*P_g_gb; 0.98*H_tes_dis; -(H_tes_cha)]'; h = bar(H1, 'stacked'); hold on; plot(1:24, H_load_h, '-r*'); grid on; colors_h = [ 0.7 0.2 0.7 % CHP heat 0.9 0.4 0.2 % Thermal heat 0.3 0.6 1 % GB heat 0.5 0.8 0.2 % TES discharge 0.8 0.3 0.3 % TES charge ]; for i = 1:size(colors_h,1) h(i).FaceColor = colors_h(i,:); h(i).EdgeColor = 'none'; end hold on; legend('CHP','CTP','GB','TES(dis)','TES(cha)','Thermal Load'); xlabel('Time (h)'); ylabel('Power (MW)'); axis([0,25,-50,290]); xticks(0:2:24); title('Thermal Power Balance'); % Nature Gas Balance subplot(1,3,3); V_chp1 = (P_g_chp(1,:)); V_gb = (P_g_gb); V_ch4 = (((P_e_p2g*0.85-P_h2_chp)*0.75)); V_buy_g = V_chp1 + V_gb - V_ch4; V = [V_chp1; V_gb; -V_ch4]'; g = bar(V, 'stacked'); hold on; plot(1:24, value(V_buy_g), '-b*'); grid on; colors_g = [ 0.7 0.2 0.7 % CHP 0.3 0.6 1 % GB 0.8 0.3 0.3 % P2G ]; hold on; for i = 1:size(colors_g,1) g(i).FaceColor = colors_g(i,:); g(i).EdgeColor = 'none'; end legend('CHP','GB','P2G','Gas Purchase'); xlabel('Time (h)'); ylabel('Gas Flow (MW)'); axis([0,25,-50,200]); xticks(0:2:24); title('Nature Gas Balance'); figure(3); E_SOC_elec(1)=SOC_ESS_max*0.3; H_SOC_heat(1)=SOC_HSS_max*0.3; for i=2:25 E_SOC_elec(i)=0.3*SOC_ESS_max+sum(P_ess_dis(1,1:i-1))/0.95+0.95*sum(P_ess_cha(1,1:i-1)); H_SOC_heat(i)=0.3*SOC_HSS_max+sum(0.98*H_tes_cha(1,1:i-1)+k_eb*P_e_eb(1,1:i-1)-H_tes_dis(1,1:i-1)); end % 子图1:电储能状态变化 subplot(2,1,1); P1 = [P_ess_cha; P_ess_dis]'; bar(P1, 'stacked'); hold on plot(0:24, E_SOC_elec, 'g*-'); xlabel('Time (h)'); ylabel('ESS (MWh)'); title('Electric Energy Storage State Variation'); grid on; ylim([-30, 60]); xticks(0:2:24); % 子图2:热储能状态变化 subplot(2,1,2); H1 = [-H_tes_cha; H_tes_dis]'; bar(H1, 'stacked'); hold on plot(0:24, H_SOC_heat, 'r*-'); xlabel('Time (h)'); ylabel('TES (MWh)'); title('Thermal Energy Storage State Variation'); grid on; ylim([-30, SOC_HSS_max]); xticks(0:2:24); % 火电机组运行轨迹 figure(4); % 定义四边形可行域顶点坐标 x_vertex = [0, 114.54, 196.36, 0]; y_vertex = [108, 86.73, 220.91, 270]; % 填充可行域(SCI素雅配色) fill(x_vertex, y_vertex,[0.8 0.9 0.8], 'EdgeColor', 'none'); hold on; % 绘制运行点轨迹(SCI专业配色,优化标记尺寸) plot(value(P_th_h), value(P_th_e), 'o-', 'Color', [0.80, 0.35, 0.10], 'MarkerSize', 6); xlabel('Heat Output (MW)'); ylabel('Power Output (MW)'); title('Thermal Power Unit Operation Trajectory'); legend('Feasible Region', 'Operation Points'); grid on; axis([-10, 210, 70, 280]); figure(5) for i=2:25 m_th_so2(i-1)=(P_th_e(1,i-1)+P_th_h(1,i-1))*6.48; m_th_nox(i-1)=(P_th_e(1,i-1)+P_th_h(1,i-1))*2.88; m_coap_so2(i-1)=-B_coap_x1(1,i-1)+m_th_so2(i-1); m_coap_nox(i-1)=-B_coap_x1(2,i-1)+m_th_nox(i-1); end level_cool=zeros(24,1); color_so2_gen = [0.80, 0.35, 0.10]; color_so2_ads = [0.90, 0.55, 0.20]; color_nox_gen = [0.10, 0.30, 0.60]; color_nox_ads = [0.20, 0.50, 0.80]; color_level = [0.8 0.9 0.8]; color_power = [0.10, 0.50, 0.20]; subplot(2,1,1); yyaxis left; plot(1:24, value(m_th_so2), '-', 'Color', color_so2_gen, 'DisplayName', 'SO₂ Generation Amount'); hold on; plot(1:24, value(m_coap_so2), '--', 'Color', color_so2_ads, 'DisplayName', 'SO₂ Adsorption Capacity'); ylabel('SO₂ mass (t)'); yyaxis right; plot(1:24, value(m_th_nox), '-', 'Color', color_nox_gen, 'DisplayName', 'NOₓ Generation Amount'); plot(1:24, value(m_coap_nox), '--', 'Color', color_nox_ads, 'DisplayName', 'NOₓ Adsorption Capacity'); ylabel('NOₓ mass (t)'); title('Pollutant Gas Treatment'); xlabel('Time (h)'); legend('Location', 'best'); grid on; xticks(0:2:24); subplot(2,1,2); yyaxis left; bar(1:24, level_cool, 0.4, 'FaceColor', color_level, 'EdgeColor', 'none', 'DisplayName', 'Operation Level (Not activated)'); ylabel('Cooling Stage'); ylim([0, 5]); yyaxis right; plot(1:24, value(P_e_coap), '-*', 'Color', color_power, 'MarkerSize', 6, 'DisplayName', 'COAP Power Consumption'); ylabel('Power (MW)'); title('COAP Operation Status'); xlabel('Time (h)'); legend('Location', 'best'); grid on; xticks(0:2:24); %% 场景2:考虑碳交易+P2G+CCS %% 初始化 close all clear clc warning off load date.mat %加载电负荷、热负荷 c=[0.0004818,2.695,52,0.0004751,0.3536,0.001004]*0.1471; p_coal=800; p_gas=3.2; %煤耗特性热电联产/火电(注:除1000,将千瓦时转化为兆瓦时)/一吨标煤价格 a=1000*[0.17,0.17,0.17,0.17,0.17,0.17,0.17,0.49,0.83,0.83,0.83,0.83,0.49,0.49,0.49,0.49,0.49,0.83,0.83,0.83,0.83,0.49,0.49,0.49]; V1=5*[10,10.4,11.3,11.7,12.4,13.2,14.1,14.6,15.1,15.3,14.7,14.4,14,13.8,13.3,12.7,12.2,12.3,11.9,11.8,11.4,11.1,10.4,10.2]; b=180; %售电售热价格 pi(1)=1; % 气象数据(辐照度 kW/m&sup2;,风速 m/s) P_Solar = [0,0,0,0,0,0.09,1.15,2.82,4.54,6.20,7.44,8.01,7.93,7.09,5.93,4.57,2.97,1.27,0.15,0,0,0,0,0]; P_Wind = [7.14,7.20,7.22,6.89,7.18,7.07,6.43,5.97,6.33,6.56,6.63,6.70,6.65,6.92,6.86,6.85,6.89,7.01,7.02,6.97,7.05,7.06,7.03,7.10]; % 光伏模型 (单位: MW) S_PV = 40000*3.75; % 光伏面积 (m&sup2;) k_eta_PV = 0.157; % 光电转换效率 P_PV_max = zeros(1,24); for i = 1:24 P_PV_max(i) = (S_PV * P_Solar(i) * k_eta_PV) / 1000; end % 风电模型 (单位: MW) Q_S_WT = 80000*4.8; % 单机额定容量 (kW) v_cut_in = 3; v_rated = 12; v_cut_out = 20; % 风速参数 P_WT_max = zeros(1,24); for i = 1:24 v = P_Wind(i); if v >= v_cut_in && v <= v_rated P_WT_max(i) = Q_S_WT * (v - v_cut_in) / (v_rated - v_cut_in) / 1000; elseif v > v_rated && v <= v_cut_out P_WT_max(i) = Q_S_WT / 1000; else P_WT_max(i) = 0; end end nt=1; %场景生成、削减,得到风光出力场景/生成nt个场景 %% 定义变量 H_th_h=sdpvar(1,24,'full'); %CHP机组热出力 P_th_e=sdpvar(1,24,'full'); %火电机组电出力 H_tes_dis=sdpvar(1,24,'full'); %储热罐储能热出力 H_tes_cha=sdpvar(1,24,'full'); %储热罐储热量 P_eb_e=sdpvar(1,24,'full'); %电锅炉出力 P_gb_g=sdpvar(1,24,'full'); %电锅炉出力 P_chp=sdpvar(1,24,'full'); %CHP机组出力 z_chp=binvar(1,24); z_chp2=binvar(1,24); %CHP机组热电状态 cr=binvar(1,24); %储热充放能状态,0充,1放; P_loss_e=sdpvar(1,24,'full'); %失负荷,5倍价格购电 m_ccs_st=sdpvar(1,24,'full'); %ccs解吸塔处理量 m_ccs_ab=sdpvar(1,24,'full'); %ccs吸收处理量 P_p2g_e=sdpvar(1,24,'full'); %P_g2p功率 P_chp_H2=sdpvar(1,24,'full'); % Y_h2=sdpvar(1,24,'full'); B=Y_h2/(1-Y_h2) % 燃气掺氢比 level_cool=binvar(3,24); % COAP冷却塔温度,进行cool级冷却(一级20、二级0、三级-20) coap_mt=sdpvar(1,24,'full'); % 活性炭质量 P_coap_e=sdpvar(1,24,'full'); % coap能耗 P_ess_e=sdpvar(1,24,'full'); % ES功率 P_pv_e=sdpvar(1,24,'full'); %光伏上网功率 P_wt_e=sdpvar(1,24,'full'); %风电上网功率 P_ess_cha=sdpvar(1,24,'full'); %储能充电功率 P_ess_dis=sdpvar(1,24,'full'); %储能放电功率 coap_x1=sdpvar(2,24,'full'); %coap辅助变量 P_load_e=date(1,:); H_h=date(2,:); H_l=date(3,:); H_load_h=H_h+H_l; SOC_HSS_max=80; SOC_ESS_max=40; E_SOC_elec(1)=SOC_ESS_max*0.3; H_SOC_heat(1)=SOC_HSS_max*0.3; Pebmax=10; Pgbmax=40; k_eb=0.99; k_gb=0.8; K_chp_H=0.4; K_chp_P=0.35; k_ccs_co2=1/(25*1.977*0.001); v_ccs0=6000; %% 约束条件 n=24; C=[]; for i=1:n C=[C,(0<=H_th_h(1,i))&(H_th_h(1,i)<=2*98.18)]; % 火电热功率约束 C=[C,implies(H_th_h(1,i)>=2*57.27,z_chp(1,i)==1)]; C=[C,z_chp(1,i)+z_chp2(1,i)==1]; C=[C,implies(z_chp2(1,i),(P_th_e(1,i)<=2*135-0.25*H_th_h(1,i))&(P_th_e(1,i)>=2*54-0.1857*H_th_h(1,i)))]; C=[C,implies(z_chp(1,i),(P_th_e(1,i)<=2*135-0.25*H_th_h(1,i))&(P_th_e(1,i)>=1.64*H_th_h(1,i)-2*50.56))]; % 火电热电耦合特性约束 C=[C,(2*43.37<=P_th_e(1,i))&(P_th_e(1,i)<=2*135)]; % 火电电功率约束 C=[C,(0<=P_chp(1,i))&(P_chp(1,i)<=2*50)]; % CHP机组总功率约束 C=[C,(0<=P_eb_e(1,i))&(P_eb_e(1,i)<=Pebmax)]; % 电锅炉出力约束 C=[C,P_loss_e(1,i)>=0]; % 失负荷惩罚为正,发电量不够5倍惩罚 C=[C,(0<=H_tes_cha(1,i))&(H_tes_cha(1,i)<=20)]; % 储热储能功率约束 C=[C,(0<=H_tes_dis(1,i))&(H_tes_dis(1,i))<=20]; % 储热释能功率约束 C=[C,(0<=0.3*SOC_HSS_max+sum(0.98*H_tes_cha(1,1:i)+k_eb*P_eb_e(1,1:i)-H_tes_dis(1,1:i)))&(0.3*SOC_HSS_max+sum(0.98*H_tes_cha(1,1:i)+k_eb*P_eb_e(1,1:i)-H_tes_dis(1,1:i))<=SOC_HSS_max)]; % 储热容量约束 C=[C,(P_th_e(1,i)+K_chp_P*P_chp(1,i)+P_wt_e(1,i)+P_pv_e(1,i)+P_loss_e(1,i)-P_eb_e(1,i)-P_p2g_e(1,i)-P_coap_e(1,i)+P_ess_e(1,i)-(5+0.269*m_ccs_st(1,i)))==P_load_e(i)]; % 电功率平衡约束 C=[C,implies(cr(i)==0,H_tes_dis(1,i)==0)]; C=[C,implies(cr(i)==1,(H_tes_cha(1,i)==0))]; % 储热不能同时充放热 C=[C,H_th_h(1,i)+0.98*H_tes_dis(1,i)-H_tes_cha(1,i)+K_chp_H*P_chp(1,i)+k_gb*P_gb_g(1,i)==H_load_h(i)]; % 热功率平衡约束 C=[C,(0<=P_gb_g(1,i))&(P_gb_g(1,i)<=Pgbmax)]; % 燃气锅炉功率约束 C=[C,(0<=m_ccs_st(1,i))&(m_ccs_st(1,i)<=120)]; % CCS再生塔最大碳捕集量限制 C=[C,(0.1*0.9*2.57*((0.4104*P_th_e(1,i)+0.06618*H_th_h(1,i)+7.092))<=m_ccs_ab(1,i))&(m_ccs_ab(1,i)<=0.9*2.57*((0.4104*P_th_e(1,i)+0.06618*H_th_h(1,i)+7.092)))]; % CCS吸收塔吸收量限制 C=[C,(0<=(v_ccs0+sum(k_ccs_co2*m_ccs_ab(1,1:i)-k_ccs_co2*m_ccs_st(1,1:i))))&(((v_ccs0+sum(k_ccs_co2*m_ccs_ab(1,1:i)-k_ccs_co2*m_ccs_st(1,1:i))))<=12000)]; % ccs富液罐体积约束 C=[C,(0<=(v_ccs0+sum(-k_ccs_co2*m_ccs_ab(1,1:i)+k_ccs_co2*m_ccs_st(1,1:i))))&(((v_ccs0+sum(-k_ccs_co2*m_ccs_ab(1,1:i)+k_ccs_co2*m_ccs_st(1,1:i))))<=12000)]; % ccs贫液罐体积约束 C=[C,(0<=P_p2g_e(1,i))&(P_p2g_e(1,i)<=30)]; % P2G功率约束 C=[C,(0<=P_chp_H2(1,i))&(P_chp_H2(1,i)<=P_p2g_e(1,i))]; % 掺氢能量约束 C=[C,(P_chp(1,i)+P_gb_g(1,i))*0.003539>=5*P_chp_H2(1,i)*0.01083]; % CHP机组掺氢比约束 C=[C,level_cool(1,i)+level_cool(2,i)+level_cool(3,i)==1]; C=[C,level_cool(3,i)==1]; % 冷却级别选择 C=[C,implies(level_cool(1,i)==1,P_coap_e(1,i)==0)]; C=[C,implies(level_cool(2,i)==1,P_coap_e(1,i)==0)]; C=[C,implies(level_cool(3,i)==1,P_coap_e(1,i)==0)]; % COAP能耗 C=[C,(0<=coap_mt(1,i))&(coap_mt(1,i)<=20)]; % 活性炭处理能力约束 C=[C,(-20<=P_ess_e(1,i))&(P_ess_e(1,i)<=20)]; % 电化学储能功率约束 C=[C,(-20<=P_ess_cha(1,i))&(P_ess_cha(1,i)<=0)]; C=[C,(0<=P_ess_dis(1,i))&(P_ess_dis(1,i)<=20)]; C=[C,implies(P_ess_e(1,i)>=0,(P_ess_dis(1,i)==P_ess_e(1,i))&(P_ess_cha(1,i)==0))]; C=[C,implies(P_ess_e(1,i)<=0,(P_ess_dis(1,i)==0)&(P_ess_cha(1,i)==P_ess_e(1,i)))]; C=[C,(0<=0.3*SOC_ESS_max+sum(P_ess_dis(1,1:i))/0.95+0.95*sum(P_ess_cha(1,1:i)))&(0.3*SOC_ESS_max+sum(P_ess_dis(1,1:i))/0.95+0.95*sum(P_ess_cha(1,1:i))<=SOC_ESS_max)]; % 电化学储能容量约束 C=[C,implies(level_cool(1,i)==1,coap_x1(1,i)==(P_th_e(1,i)+H_th_h(1,i))*6.48-0)]; C=[C,implies(level_cool(2,i)==1,coap_x1(1,i)==(P_th_e(1,i)+H_th_h(1,i))*6.48-0)]; C=[C,implies(level_cool(3,i)==1,coap_x1(1,i)==(P_th_e(1,i)+H_th_h(1,i))*6.48-0)]; C=[C,implies(level_cool(1,i)==1,coap_x1(2,i)==(P_th_e(1,i)+H_th_h(1,i))*2.88-0)]; C=[C,implies(level_cool(2,i)==1,coap_x1(2,i)==(P_th_e(1,i)+H_th_h(1,i))*2.88-0)]; C=[C,implies(level_cool(3,i)==1,coap_x1(2,i)==(P_th_e(1,i)+H_th_h(1,i))*2.88-0)]; %新能源出力约束 C=[C,(0<=P_pv_e(1,i))&(P_pv_e(1,i)<=P_PV_max(i))]; C=[C,(0<=P_wt_e(1,i))&(P_wt_e(1,i)<=P_WT_max(i))]; end for i=2:n C=[C,(2*-28.6<=P_chp(1,i)-P_chp(1,i-1))&(P_chp(1,i)-P_chp(1,i-1)<=2*28.6)]; % 燃气轮机爬坡约束 end C=[C,(-1<=sum(0.98*H_tes_cha(1,1:24)+k_eb*P_eb_e(1,1:24)-H_tes_dis(1,1:24))<=1)]; C=[C,(-1<=sum(P_ess_dis(1,1:24))/0.95+0.95*sum(P_ess_cha(1,1:24))<=1)]; % C=[C,(sum(P_es(1,1:24))==0)]; %% 目标函数 F=0; I_sale0=0; C00=0; punish0=0; C_cs=0; C_Ct0=0; mc0=0; Ec0=0; CB0=0; c_ccs0=0; for k=1:nt for i=1:n C_th1=p_coal*(c(3)+c(2)*P_th_e(1,i)+c(5)*H_th_h(1,i)+c(1)*P_th_e(1,i)^2+c(6)*P_th_e(1,i)*H_th_h(1,i)+c(4)*H_th_h(1,i)^2); % 火电机组燃料成本 I_ch4=p_gas*(((P_p2g_e(1,i)*0.85-P_chp_H2(1,i))*0.75)/0.01083); % P2G产气收益 C_chp1=p_gas*(P_chp(1,i)/0.01083); C_gb=p_gas*(P_gb_g(1,i)/0.01083); % CHP机组燃料成本 m_so2=level_cool(1,i)*coap_mt(1,i)*22.8761+level_cool(2,i)*coap_mt(1,i)*55.7041+level_cool(3,i)*coap_mt(1,i)*100.008; m_nox=level_cool(1,i)*coap_mt(1,i)*32.333+level_cool(2,i)*coap_mt(1,i)*116.349+level_cool(3,i)*coap_mt(1,i)*169.142; % 污染气体吸收量 QQ=P_th_e(1,i)+H_th_h(1,i); C_cb=6.4*(coap_x1(1,i))+8*(coap_x1(2,i))+1+2+0.16; CB(i)=C_cb; CB0=CB0+pi(k)*C_cb; % 环保代价 mc=2.57*((C_th1)/p_coal)+0.054*0.98*3.67*(P_chp(1,i)+P_gb_g(1,i)-P_chp_H2(1,i)); % 碳排放 Ec=0.69135*(P_th_e(1,i)+K_chp_P*P_chp(1,i))+0.3*(H_th_h(1,i)+K_chp_H*P_chp(1,i)+k_gb*P_gb_g(1,i)); % 碳配额 C_ct=205*(mc-Ec-m_ccs_ab(1,i)); % 碳交易 C_buy=5*a(i)*P_loss_e(1,i); H_eb=k_eb*P_eb_e(1,i); H_gb=k_gb*P_gb_g(1,i); C_cs=C_cs+pi(k)*(C_th1+C_gb+C_chp1-p_gas*P_chp_H2(1,i)/(0.01083)-I_ch4); %燃料成本 Pccs=5+0.269*m_ccs_st(1,i); % 碳捕集能耗 c_ccs=50*(m_ccs_st(1,i)-(0.054*0.98*3.67*((P_p2g_e(1,i)*0.85-P_chp_H2(1,i))*0.75))); c_ccs0=c_ccs0+pi(k)*c_ccs; % 碳封存成本 I_sale=a(i)*(P_th_e(1,i)+K_chp_P*P_chp(1,i)+P_wt_e(1,i)-P_p2g_e(1,i)+P_pv_e(1,i)-P_eb_e(1,i)-P_coap_e(1,i)-Pccs+P_ess_e(1,i))+b*(H_th_h(1,i)+0.95*H_tes_dis(1,i)-H_tes_cha(1,i)+K_chp_H*P_chp(1,i)+H_gb); %售热售电收益 C0=C_th1+C_gb+C_chp1-p_gas*P_chp_H2(1,i)/(0.01083)+C_cb+C_buy+C_ct+c_ccs-I_ch4; % 成本 I_sale0=I_sale0+pi(k)*I_sale; C00=C00+pi(k)*C0; C_Ct0=C_Ct0+pi(k)*C_ct; mc0=mc0+pi(k)*mc; Ec0=Ec0+pi(k)*Ec; F=F+pi(k)*(I_sale-C0); end end %%其他成本/固定收入 C_L=(4589000*((135*2)*(1+0.07+0.12))+100)/(365*60); %火力发电建设成本 C_com=(0.1536/(1-0.1536))*(C_L+C_cs); %火力发电运维成本 C_pv=(6000*1000*200)/(365*25)+12*sum(P_PV_max(1:24)); %光伏成本(建设加运维) C_w=(8000*1000*440)/(365*20)+20*sum(P_WT_max(1:24)); %风电成本(建设加运维) C_eb=(3200*1000*Pebmax)/(365*20)+2*sum(P_eb_e(1:24)); %电加热成本 C_yan=(300*1000*SOC_HSS_max)/(365*5)+1.5*sum(sum(H_tes_cha)); %储热罐成本 C_es=(2000*1000*SOC_ESS_max)/(365*5); %电化学储能成本 I_c=(330*1000*(370))/365; %煤电容量电价补偿 C_q=C_L+C_com+C_pv+C_w+C_eb+C_yan+C_es-I_c; C_zhejiu=C_L+(6000*1000*200)/(365*25)+(8000*1000*440)/(365*20)+(3200*1000*Pebmax)/(365*20)+(300*1000*SOC_HSS_max)/(365*5)+(2000*1000*SOC_ESS_max)/(365*5); C_yunwei=C_com+12*sum(P_PV_max(1:24))+20*sum(P_WT_max(1:24))+2*sum(P_eb_e(1:24))+1.5*sum(sum(H_tes_cha)); P_v_a=P_PV_max-P_pv_e; P_w_a=P_WT_max-P_wt_e; F_a=257.396*sum(P_w_a+P_v_a); %% 计算 F=F-C_q-F_a; F=-F; ops = sdpsettings('solver','cplex', 'verbose', 2); options.cplex.MaxTime =120; ops.cplex.timelimit = 30; aa=optimize(C,F,ops); %% 结果出具、绘图 value(-F); P_th_e=double(P_th_e); H_th_h=double(H_th_h); H_tes_dis=double(H_tes_dis); H_tes_cha=double(H_tes_cha); P_eb_e=double(P_eb_e); P_chp=double(P_chp); P_loss_e=double(P_loss_e); P_gb_g=double(P_gb_g); H_gb=k_gb*P_gb_g; C_Ct0=double(C_Ct0); mc0=double(mc0); Ec0=double(Ec0); m_ccs_st=double(m_ccs_st); m_ccs_ab=double(m_ccs_ab); Pccs=5+0.269*m_ccs_st; P_p2g_e=double(P_p2g_e); P_chp_H2=double(P_chp_H2); P_coap_e=double(P_coap_e); level_cool=double(level_cool); coap_mt=double(coap_mt); P_ess_e=double(P_ess_e); CB=double(CB); coap_x1=double(coap_x1); P_pv_e=double(P_pv_e); P_wt_e=double(P_wt_e); c_ccs0=double(c_ccs0); P_ess_cha=double(P_ess_cha); P_ess_dis=double(P_ess_dis); M_CO2_P2G=(0.54*0.98*3.67*((P_p2g_e*0.85-P_chp_H2)*0.75)); M_CCS_CO2=m_ccs_st; for i=2:25 E_SOC_elec(i)=0.3*SOC_ESS_max+sum(P_ess_dis(1,1:i-1))/0.95+0.95*sum(P_ess_cha(1,1:i-1)); H_SOC_heat(i)=0.3*SOC_HSS_max+sum(0.98*H_tes_cha(1,1:i-1)+k_eb*P_eb_e(1,1:i-1)-H_tes_dis(1,1:i-1)); % m_coap_so2(i-1)=cool(1,i-1)*mt(1,i-1)*22.8761+cool(2,i-1)*mt(1,i-1)*55.7041+cool(3,i-1)*mt(1,i-1)*100.008; % m_coap_nox(i-1)=cool(1,i-1)*mt(1,i-1)*32.333+cool(2,i-1)*mt(1,i-1)*116.349+cool(3,i-1)*mt(1,i-1)*169.142; m_th_so2(i-1)=(P_th_e(1,i-1)+H_th_h(1,i-1))*6.48; m_th_nox(i-1)=(P_th_e(1,i-1)+H_th_h(1,i-1))*2.88; m_coap_so2(i-1)=-coap_x1(1,i-1)+m_th_so2(i-1); m_coap_nox(i-1)=-coap_x1(2,i-1)+m_th_nox(i-1); end %% 结果可视化 % 新能源出力 figure(1); set(groot, 'DefaultAxesFontName', 'Times New Roman'); set(groot, 'DefaultTextFontName', 'Times New Roman'); set(groot, 'DefaultLineLineWidth', 1.5); x = 1:24; for i = 1:24 if P_PV_max(i)-P_pv_e(i)>0 PC1(i)=patch([x(i)-0.4, x(i)+0.4, x(i)+0.4, x(i)-0.4], [P_pv_e(i), P_pv_e(i), P_PV_max(i), P_PV_max(i)], 'b'); set(PC1(i), 'FaceColor', [0.97, 0.75, 0.40] , 'EdgeColor', 'none'); bj=i; end if P_WT_max(i)-P_wt_e(i)>0 PC4(i)=patch([x(i)-0.4, x(i)+0.4, x(i)+0.4, x(i)-0.4], [P_wt_e(i), P_wt_e(i), P_WT_max(i), P_WT_max(i)], 'b'); set(PC4(i), 'FaceColor', [0.40, 0.70, 0.95] , 'EdgeColor', 'none'); bj1=i; end end hold on PC2=plot(x,P_pv_e,'-rv'); set(PC2, 'Color', [0.90, 0.55, 0.20] ,'LineWidth',1); hold on PC3=plot(x,P_PV_max,'-ro'); set(PC3, 'Color', [0.80, 0.35, 0.10] ,'LineWidth',1); hold on PC5=plot(x,P_wt_e,'-rv'); set(PC5, 'Color', [0.20, 0.50, 0.80],'LineWidth',1); hold on PC6=plot(x,P_WT_max,'-ro'); set(PC6, 'Color', [0.10, 0.30, 0.60] ,'LineWidth',1); grid on; legend([PC1(bj),PC2,PC3,PC4(bj1),PC5,PC6], 'Curtailed PV Power','Actual PV Output','Maximum PV Output','Curtailed WT Power','Actual WT Output','Maximum WT Output'); xlabel('time (h)'); ylabel('Power (MW)'); title('Output Curve of New Energy Units'); ylim([-1 240]); % 功率平衡图绘制 figure(2); % 全局样式设置(SCI期刊标配) set(groot, 'DefaultAxesFontName', 'Times New Roman'); set(groot, 'DefaultTextFontName', 'Times New Roman'); set(groot, 'DefaultLineLineWidth', 1.5); % set(gcf, 'Color', 'white'); % 背景设为白色 set(gcf, 'Position', [10 10 1500 350]); % 调整图窗尺寸适配排版 % 子图1:电负荷平衡图 subplot(1,3,1); P_v_a = P_PV_max - P_pv_e; P_w_a = P_WT_max - P_wt_e; P = [P_chp*K_chp_P; P_th_e; P_pv_e; P_wt_e; -P_eb_e; -Pccs; -P_p2g_e; -P_coap_e; P_ess_cha; P_ess_dis]'; e = bar(P, 'stacked'); hold on; plot(1:24, P_load_e, '-b*', 'LineWidth', 1, 'MarkerSize', 6); % 优化标记尺寸 grid on; % 电负荷平衡图配色(SCI专业色系,保留原色系逻辑并优化对比度) colors_e = [ 0.7 0.2 0.7 % CHP机组 - 紫色 0.9 0.4 0.2 % 火电机组 - 橙色 1 0.8 0.2 % 光伏机组 - 黄色 0.5 0.8 0.2 % 风电机组 - 绿色 0.6 0.6 0.6 % 电锅炉 - 灰色 0.2 0.6 0.8 % CCS能耗 - 蓝色 0.8 0.3 0.3 % P2G能耗 - 红色 0.1 0.1 0.6 % COAP能耗 - 深蓝 0.2 0.8 0.8 % 储能充电 - 青色 0.3 0.7 0.9 % 储能放电 - 蓝绿色 ]; for i = 1:size(colors_e,1) e(i).FaceColor = colors_e(i,:); e(i).EdgeColor = 'none'; % 去除边框更素雅 end % 图例标签(英文显示,优化图例位置) legend('CHP','CTP','PV','WT','EB','CCS','P2G','COAP','ESS Charging','ESS Discharging','Electric Load', 'Location', 'best', 'FontSize', 9); xlabel('Time Interval'); ylabel('Power (MW)'); axis([0,25,-200,650]); set(gca, 'FontSize', 10); set(gca, 'GridLineStyle', ':', 'GridColor', [0.7,0.7,0.7]); % 优化网格样式 xticks(0:2:24); % X轴刻度优化,避免拥挤 title('Electric Power Balance Diagram'); % 子图2:热负荷平衡图 subplot(1,3,2); H1 = [K_chp_H*P_chp; H_th_h; H_gb; 0.98*H_tes_dis; -(H_tes_cha)]'; h = bar(H1, 'stacked'); hold on; plot(1:24, H_load_h, '-r*', 'LineWidth', 1, 'MarkerSize', 6); grid on; % 热负荷平衡图配色 colors_h = [ 0.7 0.2 0.7 % CHP供热 - 紫色 0.9 0.4 0.2 % 火电供热 - 橙色 0.3 0.6 1 % 燃气锅炉 - 蓝色 0.5 0.8 0.2 % 储热罐放热 - 绿色 0.8 0.3 0.3 % 储热罐充热 - 红色 ]; for i = 1:size(colors_h,1) h(i).FaceColor = colors_h(i,:); h(i).EdgeColor = 'none'; end legend('CHP','CTP','GB','TES Discharge','TES Charging','Thermal Load', 'Location', 'best'); xlabel('Time Interval'); ylabel('Power (MW)'); axis([0,25,-50,290]); set(gca, 'FontSize', 10); set(gca, 'GridLineStyle', ':', 'GridColor', [0.7,0.7,0.7]); xticks(0:2:24); title('Thermal Power Balance Diagram'); % 子图3:气平衡图 subplot(1,3,3); V_chp1 = (P_chp(1,:)); V_gb = (P_gb_g); V_ch4 = (((P_p2g_e*0.85-P_chp_H2)*0.75)); V_buy = V_chp1 + V1 + V_gb - V_ch4; V = [-V_chp1; -V_gb; V_ch4; V_buy]'; g = bar(V, 'stacked'); hold on; plot(1:24, V1, '-m*', 'MarkerSize', 6); grid on; % 气平衡图配色 colors_g = [ 0.7 0.2 0.7 % CHP - 紫色 0.3 0.6 1 % 燃气锅炉 - 蓝色 0.8 0.3 0.3 % P2G - 红色 0.5 0.8 0.2 % 气负荷 - 绿色 ]; for i = 1:size(colors_g,1) g(i).FaceColor = colors_g(i,:); g(i).EdgeColor = 'none'; end legend('CHP','GB','P2G','Gas Purchase','Gas Load', 'Location', 'best', 'FontSize', 9); xlabel('Time Interval'); ylabel('Gas Volume (MW equivalent)'); % 补充气维度标注,更规范 axis([0,25,-150,270]); set(gca); set(gca, 'GridLineStyle', ':', 'GridColor', [0.7,0.7,0.7]); xticks(0:2:24); title('Gas Power Balance Diagram'); % 储能状态变化 figure(3); % 全局样式设置(SCI期刊标配) set(groot, 'DefaultAxesFontName', 'Times New Roman'); set(groot, 'DefaultTextFontName', 'Times New Roman'); set(groot, 'DefaultLineLineWidth', 1.5); % set(gcf, 'Color', 'white'); % 背景设为白色 % 子图1:电储能状态变化 subplot(2,1,1); % 绘制电储能曲线(SCI专业配色,优化标记尺寸) plot(0:24, E_SOC_elec, '*-', 'Color', [0.10, 0.30, 0.60], 'MarkerSize', 6); xlabel('Time (h)'); ylabel('Stored Electric Energy (MWh)'); title('Electric Energy Storage State Variation'); grid on; % 优化网格样式(浅灰色虚线,不喧宾夺主) set(gca, 'GridLineStyle', ':', 'GridColor', [0.7, 0.7, 0.7]); set(gca, 'LineWidth', 1); % 坐标轴线条宽度 ylim([0, SOC_ESS_max]); set(gca, 'FontSize', 10); xticks(0:2:24); % X轴刻度优化,避免拥挤 % 子图2:热储能状态变化 subplot(2,1,2); % 绘制热储能曲线(SCI专业配色,优化标记尺寸) plot(0:24, H_SOC_heat, '*-', 'Color', [0.80, 0.35, 0.10], 'MarkerSize', 6); xlabel('Time (h)'); ylabel('Stored Thermal Energy (MWh)'); title('Thermal Energy Storage State Variation'); grid on; set(gca, 'GridLineStyle', ':', 'GridColor', [0.7, 0.7, 0.7]); set(gca, 'LineWidth', 1); ylim([0, SOC_HSS_max]); set(gca, 'FontSize', 10); xticks(0:2:24); % 火电机组运行轨迹 figure(4); set(groot, 'DefaultAxesFontName', 'Times New Roman'); set(groot, 'DefaultTextFontName', 'Times New Roman'); set(groot, 'DefaultLineLineWidth', 1.5); % set(gcf, 'Color', 'white'); % 背景设为白色 % 定义四边形可行域顶点坐标 x_vertex = [0, 114.54, 196.36, 0]; y_vertex = [108, 86.73, 220.91, 270]; % 填充可行域(SCI素雅配色) fill(x_vertex, y_vertex,[0.8 0.9 0.8], 'EdgeColor', 'none'); hold on; % 绘制运行点轨迹(SCI专业配色,优化标记尺寸) plot(H_th_h, P_th_e, '*-', 'Color', [0.80, 0.35, 0.10], 'MarkerSize', 6); xlabel('Heat Output (MW)'); ylabel('Power Output (MW)'); title('Thermal Power Unit Operation Trajectory'); legend('Feasible Region', 'Operation Points', 'Location', 'best'); grid on; set(gca, 'GridLineStyle', ':', 'GridColor', [0.7, 0.7, 0.7]); set(gca, 'LineWidth', 1); axis([-10, 210, 70, 280]); % P2G-CCS系统分析 figure(5); set(groot, 'DefaultAxesFontName', 'Times New Roman'); % 统一字体为SCI标配 set(groot, 'DefaultTextFontName', 'Times New Roman'); set(groot, 'DefaultLineLineWidth', 1.5); % 线条宽度适配SCI印刷 % set(gcf, 'Color', 'white'); % 背景设为白色 subplot(2,1,1); yyaxis left; plot(1:24, P_p2g_e, 'b-*', 'MarkerSize', 6); % 绘制电功率输入曲线 ylabel('Electric Power Input (MW)'); ylim([-1 50]); yyaxis right; plot(1:24, V_ch4, 'r-o', 'MarkerSize', 6); % 绘制天然气产量曲线 ylabel('Natural Gas Production (10^4 m^3)'); ylim([-1 50]); title('P2G System Performance'); xlabel('Time (h)'); legend('Power Consumption', 'Gas Production', 'Location', 'best'); grid on; set(gca, 'GridLineStyle', ':', 'GridColor', [0.7,0.7,0.7]); % 网格样式优化 set(gca, 'LineWidth', 1); xticks(0:2:24); % X轴刻度优化 ax1 = gca; % 左侧Y轴 ax1.YAxis(1).Color = 'k'; ax1.YAxis(1).TickLabelColor = 'k'; ax1.YAxis(1).Label.Color = 'k'; % 右侧Y轴 ax1.YAxis(2).Color = 'k'; ax1.YAxis(2).TickLabelColor = 'k'; ax1.YAxis(2).Label.Color = 'k'; % 子图2:碳源利用分析 subplot(2,1,2); carbon_flow = [M_CCS_CO2;M_CO2_P2G]; bar(carbon_flow', 1.2, 'FaceColor', 'flat'); % 绘制碳流量柱状图 colormap([0.10,0.30,0.60; 0.80,0.35,0.10]); title('Carbon Source Utilization Analysis'); xlabel('Time (h)'); ylabel('CO_2 (ton)'); legend('CCS Capture Amount', 'P2G Consumption Amount', 'Location', 'best'); ylim([-1 170]); grid on; set(gca, 'GridLineStyle', ':', 'GridColor', [0.7,0.7,0.7]); set(gca, 'LineWidth', 1); xticks(0:2:24); % COAP处理污染物系统分析 figure(6) set(groot, 'DefaultAxesFontName', 'Times New Roman'); % 统一字体为SCI标配 set(groot, 'DefaultTextFontName', 'Times New Roman'); set(groot, 'DefaultLineLineWidth', 1.5); % 线条宽度适配SCI印刷 % set(gcf, 'Color', 'white'); % 背景设为白色 level_cool=zeros(24,1); color_so2_gen = [0.80, 0.35, 0.10]; color_so2_ads = [0.90, 0.55, 0.20]; color_nox_gen = [0.10, 0.30, 0.60]; color_nox_ads = [0.20, 0.50, 0.80]; color_level = [0.85, 0.85, 0.85]; color_power = [0.10, 0.50, 0.20]; subplot(2,1,1); yyaxis left; plot(1:24, m_th_so2, '-', 'Color', color_so2_gen, 'DisplayName', 'SO₂ Generation Amount'); hold on; plot(1:24, m_coap_so2, '--', 'Color', color_so2_ads, 'DisplayName', 'SO₂ Adsorption Capacity'); ylabel('SO₂ mass (t)'); ylim([0, max([m_th_so2, m_coap_so2,m_th_nox, m_coap_nox])*1.2]); yyaxis right; plot(1:24, m_th_nox, '-', 'Color', color_nox_gen, 'DisplayName', 'NOₓ Generation Amount'); plot(1:24, m_coap_nox, '--', 'Color', color_nox_ads, 'DisplayName', 'NOₓ Adsorption Capacity'); ylabel('NOₓ mass (t)'); ylim([0, max([m_th_so2, m_coap_so2,m_th_nox, m_coap_nox])*1.2]); title('Pollutant Gas Treatment'); xlabel('Time (h)'); legend('Location', 'best'); grid on; set(gca, 'GridLineStyle', ':', 'GridColor', [0.7,0.7,0.7]); set(gca, 'LineWidth', 1); xticks(0:2:24); % 设置子图1纵坐标刻度和字体为黑色 ax1 = gca; ax1.YAxis(1).Color = 'k'; ax1.YAxis(1).TickLabelColor = 'k'; ax1.YAxis(1).Label.Color = 'k'; ax1.YAxis(2).Color = 'k'; ax1.YAxis(2).TickLabelColor = 'k'; ax1.YAxis(2).Label.Color = 'k'; subplot(2,1,2); yyaxis left; bar(1:24, level_cool, 0.4, 'FaceColor', color_level, 'EdgeColor', 'none', 'DisplayName', 'Operation Level'); ylabel('Cooling Stage'); ylim([0, 5]); yyaxis right; plot(1:24, P_coap_e, '-*', 'Color', color_power, 'MarkerSize', 6, 'DisplayName', 'COAP Power Consumption'); ylabel('Power (MW)'); ylim([0, max(P_coap_e+100)*1.3]); title('COAP Operation Status'); xlabel('Time (h)'); legend('Location', 'best'); grid on; set(gca, 'GridLineStyle', ':', 'GridColor', [0.7,0.7,0.7]); set(gca, 'LineWidth', 1); xticks(0:2:24); % 设置子图2纵坐标刻度和字体为黑色 ax2 = gca; ax2.YAxis(1).Color = 'k'; ax2.YAxis(1).TickLabelColor = 'k'; ax2.YAxis(1).Label.Color = 'k'; ax2.YAxis(2).Color = 'k'; ax2.YAxis(2).TickLabelColor = 'k'; ax2.YAxis(2).Label.Color = 'k'; C00=double(C00); F1=[double(I_sale0)-double(punish0)+sum(V1*p_gas/0.01083),I_c,-double(C_cs)-sum(V1*p_gas/0.01083),-double(C_zhejiu),-double(C_yunwei)-c_ccs0,-C_Ct0,-sum(double(CB0)),-double(F_a),-5*sum(P_loss_e.*a)]; sum(F1) disp([' 总收益: ', num2str(sum(F1)/(1000*7.06), '%.2f'), ' 千美元']); disp('========== 运行收益分解 =========='); disp([' 售能收益: ', num2str(F1(1)/(1000*7.06), '%.2f'), ' 千美元']); disp([' 容量电价补偿收益: ', num2str(F1(2)/(1000*7.06), '%.2f'), ' 千美元']); disp([' 燃料成本: ', num2str(F1(3)/(1000*7.06), '%.2f'), ' 千美元']); disp([' 设备折旧成本: ', num2str(F1(4)/(1000*7.06), '%.2f'), ' 千美元']); disp([' 运维成本: ', num2str(F1(5)/(1000*7.06), '%.2f'), ' 千美元']); disp([' 碳交易收益/成本: ', num2str(F1(6)/(1000*7.06), '%.2f'), ' 千美元']); disp([' 环保代价处理成本: ', num2str(F1(7)/(1000*7.06), '%.2f'), ' 千美元']); disp([' 弃电惩罚成本: ', num2str(F1(8)/(1000*7.06), '%.2f'), ' 千美元']); disp([' 失负荷惩罚成本: ', num2str(F1(9)/(1000*7.06), '%.2f'), ' 千美元']); disp([' 弃风率: ',num2str(100*(1-sum(P_wt_e)/sum(P_WT_max)), '%.2f'),'%']); disp([' 弃光率: ',num2str(100*(1-sum(P_pv_e)/sum(P_PV_max)), '%.2f'),'%']); disp([' 失负荷电量: ',num2str(sum(P_loss_e), '%.2f'),'MW']); % 成本收益可视化 figure(7); conv = 1/(1000*7.06); % 转换系数:千美元 threshold = 1e-3; % 近0阈值,绝对值小于该值的项目过滤 item_names = { 'Energy Sales Revenue'; % 1:售能收益 'Capacity Tariff Revenue'; % 2:容量电价补偿收益 'Fuel Cost'; % 3:燃料成本 'Equipment Depreciation Cost'; % 4:设备折旧成本 'O&M Cost'; % 5:运维成本 'Carbon Trading Revenue/Cost'; % 6:碳交易收益/成本 'Environmental Treatment Cost'; % 7:环保代价处理成本 'Curtailed Power Penalty Cost'; % 8:弃电惩罚成本 'Load Loss Penalty Cost' % 9:失负荷惩罚成本 }; item_values = F1 * conv; % 转换为千美元 % 收益项目(正数值,过滤近0) revenue_mask = (item_values > 0) & (abs(item_values) > threshold); revenue_names = item_names(revenue_mask); revenue_values = item_values(revenue_mask); % 成本项目(负数值,过滤近0,取绝对值用于饼图) cost_mask = (item_values < 0) & (abs(item_values) > threshold); cost_names = item_names(cost_mask); cost_values = abs(item_values(cost_mask)); set(groot, 'DefaultAxesFontName', 'Times New Roman'); set(groot, 'DefaultTextFontName', 'Times New Roman'); set(groot, 'DefaultLineLineWidth', 1.5); pie_colors = [0.80,0.35,0.10; 0.10,0.30,0.60; 0.10,0.50,0.20; 0.70,0.20,0.70; 0.60,0.60,0.60; 0.20,0.60,0.80; 0.30,0.70,0.90]; bar_colors = [0.10,0.30,0.60; 0.80,0.35,0.10; 0.10,0.50,0.20]; subplot(1,2,1); explode = zeros(1, length(cost_values)); min_sector_idx = find(cost_values / sum(cost_values) < 0.02); if ~isempty(min_sector_idx) explode(min_sector_idx) = 0.5; % 极小扇区爆炸 else explode(end) = 1; end h_pie = pie(cost_values, explode); pie_patches = findobj(h_pie, 'Type', 'patch'); % 获取饼图扇区对象 for i = 1:length(pie_patches) pie_patches(i).FaceColor = pie_colors(i,:); end legend(cost_names, 'Location', 'best', 'FontSize', 8); title('Cost Breakdown (Thousand USD)', 'FontSize', 12); set(gca, 'FontSize', 10); t1 = title('Cost Breakdown (Thousand USD)', 'FontSize', 12); pos = t1.Position; t1.Position = [pos(1), pos(2) + 0.15, pos(3)]; % 增大Y坐标实现标题上移 subplot(1,2,2); h = barh(revenue_values); set(gca, 'YTick', 1:length(revenue_names)); set(gca, 'YTickLabel', revenue_names); title('Revenue Breakdown (Thousand USD)', 'FontSize', 12); xlabel('Revenue (Thousand USD)', 'FontSize', 10); grid on; set(gca, 'FontSize', 10); for i = 1:length(revenue_values) text(revenue_values(i), i, sprintf('%.2f', revenue_values(i)), ... 'HorizontalAlignment', 'left', 'VerticalAlignment', 'middle', ... 'FontSize', 9, 'FontName', 'Times New Roman'); end set(gcf, 'Position', [100 100 1200 500]); % 调整图窗大小 set(gca, 'FontSize', 10);
最新发布
12-25
评论
成就一亿技术人!
拼手气红包6.0元
还能输入1000个字符
 
红包 添加红包
表情包 插入表情
 条评论被折叠 查看
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

Claude的羽毛

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值